﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace g4
{
    /// <summary>
    /// 网格等值线计算类：在网格上生成函数值等于特定阈值（默认0）的等值线
    /// 核心原理：通过判断网格三角形与阈值的相交情况，计算交点并连接成连续曲线
    /// </summary>
    public class MeshIsoCurves
    {
        /// <summary>
        /// 目标网格（用于计算等值线的基础网格）
        /// </summary>
        public DMesh3 Mesh;

        /// <summary>
        /// 计算顶点函数值的函数（输入顶点坐标，输出对应的函数值）
        /// 等值线为该函数值等于0的点集
        /// </summary>
        public Func<Vector3d, double> ValueF = null;

        /// <summary>
        /// 可选：通过顶点ID获取函数值的函数（可使用预计算值提高效率）
        /// </summary>
        public Func<int, double> VertexValueF = null;

        /// <summary>
        /// 若为true，将内部多线程预计算所有顶点的函数值（存储于VertexValueF）
        /// </summary>
        public bool PrecomputeVertexValues = false;


        /// <summary>
        /// 求根算法模式：用于计算边与等值线的交点
        /// </summary>
        public enum RootfindingModes
        {
            SingleLerp = 0,   // 单次线性插值（快速但精度较低）
            LerpSteps = 1,    // 多步线性插值（迭代优化精度）
            Bisection = 2     // 二分法（精度高但计算成本高）
        }

        /// <summary>
        /// 选择用于计算边与等值线交点的求根方法
        /// </summary>
        public RootfindingModes RootMode = RootfindingModes.SingleLerp;

        /// <summary>
        /// 求根算法的迭代步数（SingleLerp模式忽略此参数）
        /// </summary>
        public int RootModeSteps = 5;


        /// <summary>
        /// 存储计算得到的等值线图（顶点为等值线上的点，边为连接的线段）
        /// </summary>
        public DGraph3 Graph = null;

        /// <summary>
        /// 三角形与等值线相交的类型枚举
        /// </summary>
        public enum TriangleCase
        {
            EdgeEdge = 1,   // 等值线穿过三角形的两条边（边-边相交）
            EdgeVertex = 2, // 等值线穿过三角形的一条边和一个顶点（边-顶点相交）
            OnEdge = 3      // 等值线沿三角形的一条边（边在等值线上）
        }

        /// <summary>
        /// 是否需要记录图中边的详细信息（如对应的网格三角形、边等）
        /// </summary>
        public bool WantGraphEdgeInfo = false;

        /// <summary>
        /// 图中边的详细信息结构
        /// </summary>
        public struct GraphEdgeInfo
        {
            public TriangleCase caseType;   // 相交类型
            public int mesh_tri;            // 对应的网格三角形ID
            public Index2i mesh_edges;      // 对应的网格边/顶点（根据类型不同：EdgeEdge为两条边，EdgeVertex为边和顶点，OnEdge为边和-1）
            public Index2i order;           // 边的顶点顺序（用于方向判断）
        }

        /// <summary>
        /// 图中所有边的详细信息（与Graph.Edges对应）
        /// </summary>
        public DVector<GraphEdgeInfo> GraphEdges = null;

        /// <summary>
        /// 记录边与等值线的交点位置（键为网格边ID，值为交点坐标）
        /// </summary>
        Dictionary<int, Vector3d> EdgeLocations = new Dictionary<int, Vector3d>();


        /// <summary>
        /// 构造函数：初始化网格等值线计算器
        /// </summary>
        /// <param name="mesh">目标网格</param>
        /// <param name="valueF">计算顶点函数值的函数（用于确定等值线）</param>
        public MeshIsoCurves(DMesh3 mesh, Func<Vector3d, double> valueF)
        {
            Mesh = mesh;
            ValueF = valueF;
        }

        /// <summary>
        /// 计算整个网格的等值线
        /// </summary>
        public void Compute()
        {
            compute_full(Mesh.TriangleIndices(), true);
        }

        /// <summary>
        /// 计算指定三角形集合的等值线
        /// </summary>
        /// <param name="Triangles">待处理的三角形ID集合</param>
        public void Compute(IEnumerable<int> Triangles)
        {
            compute_full(Triangles);
        }


        /* 内部实现 */

        /// <summary>
        /// 存储图中顶点坐标与ID的映射（避免重复添加顶点）
        /// </summary>
        Dictionary<Vector3d, int> Vertices;

        /// <summary>
        /// 核心计算方法：遍历三角形，判断与等值线的相交情况并构建等值线图
        /// </summary>
        /// <param name="Triangles">待处理的三角形ID集合</param>
        /// <param name="bIsFullMeshHint">是否为整个网格（用于优化预处理）</param>
        protected void compute_full(IEnumerable<int> Triangles, bool bIsFullMeshHint = false)
        {
            // 初始化等值线图和边信息存储
            Graph = new DGraph3();
            if (WantGraphEdgeInfo)
                GraphEdges = new DVector<GraphEdgeInfo>();

            // 初始化顶点映射（确保坐标相同的点只添加一次）
            Vertices = new Dictionary<Vector3d, int>();


            // 预计算顶点函数值（多线程加速）
            double[] vertexValues = null;
            if (PrecomputeVertexValues)
            {
                vertexValues = new double[Mesh.MaxVertexID];
                IEnumerable<int> verts = Mesh.VertexIndices();
                // 若只处理部分三角形，仅预计算这些三角形的顶点值
                if (!bIsFullMeshHint)
                {
                    MeshVertexSelection vertices = new MeshVertexSelection(Mesh);
                    vertices.SelectTriangleVertices(Triangles);
                    verts = vertices;
                }
                // 并行计算每个顶点的函数值
                gParallel.ForEach(verts, (vid) => {
                    vertexValues[vid] = ValueF(Mesh.GetVertex(vid));
                });
                // 将预计算值关联到VertexValueF
                VertexValueF = (vid) => { return vertexValues[vid]; };
            }


            // 遍历每个三角形，计算与等值线的交点并构建曲线
            foreach (int tid in Triangles)
            {
                // 获取三角形的三个顶点坐标和ID
                Vector3dTuple3 tv = new Vector3dTuple3();
                Mesh.GetTriVertices(tid, ref tv.V0, ref tv.V1, ref tv.V2);
                Index3i triVerts = Mesh.GetTriangle(tid);

                // 获取三个顶点的函数值（使用预计算值或实时计算）
                Vector3d f = (VertexValueF != null) ?
                    new Vector3d(VertexValueF(triVerts.a), VertexValueF(triVerts.b), VertexValueF(triVerts.c))
                    : new Vector3d(ValueF(tv.V0), ValueF(tv.V1), ValueF(tv.V2));

                // 若三角形三个顶点的函数值全为正或全为负，与等值线（0）无交点，跳过
                if (f.x < 0 && f.y < 0 && f.z < 0)
                    continue;
                if (f.x > 0 && f.y > 0 && f.z > 0)
                    continue;

                // 获取三角形的三条边ID
                Index3i triEdges = Mesh.GetTriEdges(tid);

                // 处理三角形有顶点在等值线上（函数值为0）的情况
                if (f.x * f.y * f.z == 0)
                {
                    // 找到函数值为0的顶点（z0）
                    int z0 = (f.x == 0) ? 0 : ((f.y == 0) ? 1 : 2);
                    int i1 = (z0 + 1) % 3, i2 = (z0 + 2) % 3;

                    // 若另外两个顶点的函数值同号（都正或都负），则无交点（仅顶点在等值线上），跳过
                    if (f[i1] * f[i2] > 0)
                        continue;

                    // 情况1：有两个顶点在等值线上（边在等值线上）
                    if (f[i1] == 0 || f[i2] == 0)
                    {
                        int z1 = f[i1] == 0 ? i1 : i2;
                        // 调整顶点顺序以保持一致性
                        if ((z0 + 1) % 3 != z1)
                        {
                            int tmp = z0; z0 = z1; z1 = tmp;
                        }
                        // 添加两个顶点到图中并连接
                        int e0 = add_or_append_vertex(Mesh.GetVertex(triVerts[z0]));
                        int e1 = add_or_append_vertex(Mesh.GetVertex(triVerts[z1]));
                        int graphEid = Graph.AppendEdge(e0, e1, (int)TriangleCase.OnEdge);
                        // 记录边信息（若需要）
                        if (graphEid >= 0 && WantGraphEdgeInfo)
                            add_on_edge(graphEid, tid, triEdges[z0], new Index2i(e0, e1));
                    }
                    // 情况2：只有一个顶点在等值线上，另一条边穿过等值线（边-顶点相交）
                    else
                    {
                        Util.gDevAssert(f[i1] * f[i2] < 0); // 确保另外两个顶点异号

                        // 添加等值线上的顶点
                        int vertVid = add_or_append_vertex(Mesh.GetVertex(triVerts[z0]));

                        // 计算边（i1-i2）与等值线的交点
                        int i = i1, j = i2;
                        if (triVerts[j] < triVerts[i])
                        {
                            int tmp = i; i = j; j = tmp; // 保持顶点顺序一致
                        }
                        Vector3d cross = find_crossing(tv[i], tv[j], f[i], f[j]);
                        int crossVid = add_or_append_vertex(cross);
                        add_edge_pos(triVerts[i], triVerts[j], cross); // 记录边的交点

                        // 连接顶点和交点（若不重合）
                        if (vertVid != crossVid)
                        {
                            int graphEid = Graph.AppendEdge(vertVid, crossVid, (int)TriangleCase.EdgeVertex);
                            // 记录边信息（若需要）
                            if (graphEid >= 0 && WantGraphEdgeInfo)
                                add_edge_vert(graphEid, tid, triEdges[(z0 + 1) % 3], triVerts[z0], new Index2i(vertVid, crossVid));
                        }
                    }
                }
                // 处理三角形无顶点在等值线上，但有边穿过等值线（边-边相交）
                else
                {
                    Index3i crossVerts = Index3i.Min; // 存储交点对应的图顶点ID
                    int lessThan = 0; // 函数值小于0的顶点数量

                    // 遍历三角形的三条边，计算与等值线的交点
                    for (int tei = 0; tei < 3; ++tei)
                    {
                        int i = tei, j = (tei + 1) % 3;
                        if (f[i] < 0) lessThan++; // 统计负函数值顶点数量
                        if (f[i] * f[j] > 0) continue; // 边的两个顶点同号，无交点

                        // 调整顶点顺序以保持一致性
                        if (triVerts[j] < triVerts[i])
                        {
                            int tmp = i; i = j; j = tmp;
                        }
                        // 计算交点并添加到图中
                        Vector3d cross = find_crossing(tv[i], tv[j], f[i], f[j]);
                        crossVerts[tei] = add_or_append_vertex(cross);
                        add_edge_pos(triVerts[i], triVerts[j], cross); // 记录边的交点
                    }

                    // 确定有效交点的索引（等值线穿过两条边，对应两个交点）
                    int e0 = (crossVerts.a == int.MinValue) ? 1 : 0;
                    int e1 = (crossVerts.c == int.MinValue) ? 1 : 2;
                    if (e0 == 0 && e1 == 2)
                    {
                        e0 = 2; e1 = 0; // 保持方向一致性
                    }

                    // 根据负函数值顶点数量调整方向（确保等值线方向一致）
                    if (lessThan == 1)
                    {
                        int tmp = e0; e0 = e1; e1 = tmp;
                    }

                    // 连接两个交点（若不重合）
                    int ev0 = crossVerts[e0];
                    int ev1 = crossVerts[e1];
                    if (ev0 != ev1)
                    {
                        Util.gDevAssert(ev0 != int.MinValue && ev1 != int.MinValue);
                        int graphEid = Graph.AppendEdge(ev0, ev1, (int)TriangleCase.EdgeEdge);
                        // 记录边信息（若需要）
                        if (graphEid >= 0 && WantGraphEdgeInfo)
                            add_edge_edge(graphEid, tid, new Index2i(triEdges[e0], triEdges[e1]), new Index2i(ev0, ev1));
                    }
                }
            }

            // 释放顶点映射（避免占用内存）
            Vertices = null;
        }


        /// <summary>
        /// 向图中添加顶点（若已存在则返回现有ID，否则新建）
        /// </summary>
        /// <param name="pos">顶点坐标</param>
        /// <returns>图中顶点ID</returns>
        int add_or_append_vertex(Vector3d pos)
        {
            int vid;
            if (!Vertices.TryGetValue(pos, out vid))
            {
                vid = Graph.AppendVertex(pos);
                Vertices.Add(pos, vid);
            }
            return vid;
        }


        /// <summary>
        /// 记录OnEdge类型的边信息
        /// </summary>
        void add_on_edge(int graphEid, int meshTri, int meshEdge, Index2i order)
        {
            GraphEdgeInfo einfo = new GraphEdgeInfo()
            {
                caseType = TriangleCase.OnEdge,
                mesh_edges = new Index2i(meshEdge, -1),
                mesh_tri = meshTri,
                order = order
            };
            GraphEdges.insertAt(einfo, graphEid);
        }

        /// <summary>
        /// 记录EdgeVertex类型的边信息
        /// </summary>
        void add_edge_vert(int graphEid, int meshTri, int meshEdge, int meshVert, Index2i order)
        {
            GraphEdgeInfo einfo = new GraphEdgeInfo()
            {
                caseType = TriangleCase.EdgeVertex,
                mesh_edges = new Index2i(meshEdge, meshVert),
                mesh_tri = meshTri,
                order = order
            };
            GraphEdges.insertAt(einfo, graphEid);
        }

        /// <summary>
        /// 记录EdgeEdge类型的边信息
        /// </summary>
        void add_edge_edge(int graphEid, int meshTri, Index2i meshEdges, Index2i order)
        {
            GraphEdgeInfo einfo = new GraphEdgeInfo()
            {
                caseType = TriangleCase.EdgeEdge,
                mesh_edges = meshEdges,
                mesh_tri = meshTri,
                order = order
            };
            GraphEdges.insertAt(einfo, graphEid);
        }


        /// <summary>
        /// 计算边（a-b）与等值线（函数值0）的交点
        /// </summary>
        /// <param name="a">边的起点坐标</param>
        /// <param name="b">边的终点坐标</param>
        /// <param name="fA">a点的函数值</param>
        /// <param name="fB">b点的函数值</param>
        /// <returns>交点坐标</returns>
        Vector3d find_crossing(Vector3d a, Vector3d b, double fA, double fB)
        {
            // 确保fA <= 0且fB >= 0（调整点的顺序）
            if (fB < fA)
            {
                Vector3d tmp = a; a = b; b = tmp;
                double f = fA; fA = fB; fB = f;
            }

            // 二分法求根（高精度）
            if (RootMode == RootfindingModes.Bisection)
            {
                for (int k = 0; k < RootModeSteps; ++k)
                {
                    Vector3d c = Vector3d.Lerp(a, b, 0.5); // 中点
                    double f = ValueF(c);
                    if (f < 0)
                    {
                        fA = f; a = c; // 交点在c-b之间
                    }
                    else
                    {
                        fB = f; b = c; // 交点在a-c之间
                    }
                }
                return Vector3d.Lerp(a, b, 0.5); // 返回最终区间的中点
            }
            // 线性插值法（快速）
            else
            {
                // 避免除以零（函数值无变化，返回起点）
                if (Math.Abs(fB - fA) < MathUtil.ZeroTolerance)
                    return a;

                double t = 0;
                // 多步线性插值（迭代优化精度）
                if (RootMode == RootfindingModes.LerpSteps)
                {
                    for (int k = 0; k < RootModeSteps; ++k)
                    {
                        // 线性插值计算t（f(t) = 0时的参数）
                        t = MathUtil.Clamp((0 - fA) / (fB - fA), 0, 1);
                        Vector3d c = (1 - t) * a + t * b;
                        double f = ValueF(c);
                        if (f < 0)
                        {
                            fA = f; a = c; // 交点在c-b之间
                        }
                        else
                        {
                            fB = f; b = c; // 交点在a-c之间
                        }
                    }
                }

                // 计算最终插值参数并返回交点
                t = MathUtil.Clamp((0 - fA) / (fB - fA), 0, 1);
                return (1 - t) * a + t * b;
            }
        }


        /// <summary>
        /// 记录边与等值线的交点位置（避免重复记录）
        /// </summary>
        /// <param name="a">边的起点ID</param>
        /// <param name="b">边的终点ID</param>
        /// <param name="crossingPos">交点坐标</param>
        void add_edge_pos(int a, int b, Vector3d crossingPos)
        {
            int eid = Mesh.FindEdge(a, b);
            if (eid == DMesh3.InvalidID)
                throw new Exception("MeshIsoCurves.add_edge_split: 无效的边ID");
            if (!EdgeLocations.ContainsKey(eid))
                EdgeLocations[eid] = crossingPos;
        }


        /// <summary>
        /// 在网格上沿等值线交点分割边（忽略过短的边或过近的顶点）
        /// </summary>
        /// <param name="minLen">最小边长度（小于此值的边不分割）</param>
        public void SplitAtIsoCrossings(double minLen = 0)
        {
            foreach (var pair in EdgeLocations)
            {
                int eid = pair.Key;
                Vector3d pos = pair.Value;
                if (!Mesh.IsEdge(eid))
                    continue;

                // 检查边的长度和交点位置（避免过短边或过近顶点）
                Index2i ev = Mesh.GetEdgeV(eid);
                Vector3d a = Mesh.GetVertex(ev.a);
                Vector3d b = Mesh.GetVertex(ev.b);
                if (a.Distance(b) < minLen)
                    continue;
                Vector3d mid = (a + b) * 0.5;
                if (a.Distance(mid) < minLen || b.Distance(mid) < minLen)
                    continue;

                // 分割边并设置新顶点位置为交点
                DMesh3.EdgeSplitInfo splitInfo;
                if (Mesh.SplitEdge(eid, out splitInfo) == MeshResult.Ok)
                    Mesh.SetVertex(splitInfo.vNew, pos);
            }
        }


        /// <summary>
        /// 判断图中边的方向是否需要反转以与网格方向一致
        /// 主要用于DGraph3Util.ExtractCurves提取曲线时保持方向一致性
        /// </summary>
        /// <param name="graphEid">图中边的ID</param>
        /// <returns>若需要反转则返回true，否则返回false</returns>
        public bool ShouldReverseGraphEdge(int graphEid)
        {
            if (GraphEdges == null)
                throw new Exception("MeshIsoCurves.OrientEdge: 需启用WantGraphEdgeInfo以记录边信息");

            Index2i graphEv = Graph.GetEdgeV(graphEid);
            GraphEdgeInfo einfo = GraphEdges[graphEid];

            // 若图中边的顶点顺序与记录的顺序相反，则需要反转
            if (graphEv.b == einfo.order.a && graphEv.a == einfo.order.b)
            {
                return true;
            }
            Util.gDevAssert(graphEv.a == einfo.order.a && graphEv.b == einfo.order.b);
            return false;
        }
    }
}